% SCRIPT THAT GENERATES A MOVING-BOUNDARY DOMAIN WHOSE VELOCITY
% ACCELERATES/DECELERATES IN THE SECOND HALF OF THE PERIOD. 
% THE SPEED OF CONTRACTION/EXPANSION IN THE FIRST HALF OF THE PERIOD IS CONTROLLED BY om1 IN THE SECOND HALF BY  om2, 
% A PHASE phi1 phi2 IS ADDED TO THE FORCING AND MATCHED EVERY TIME THE VELOCITY SWITCHES
% SO TO GUARANTEE THE CONTINUITY OF THE COORDINATES
clear all
close all
% GRID POINTS
n=160;
% RADIOUS OF THE SPHERE
R=2.0;
% VISCOSITY
mu=1;
% FREQUENCY OF THE BOUNDARY
om=1;
om1=om;
om2=om/15;
phi1=0;
% INITIAL TIME
t=0;
% FINAL TIME
t_fin=2*pi/om;
% TIME STEP
dt0=2*pi*om/(100);
% BOUNDARY AT TIME t=0
phi=phi1;
[xb,yb,ub,vb,theta,dtheta]=moving_coord2D_phi(n,R,t,om1,phi);
% INIZIALIZATIONS
i=0;
modom_old=0;

while (om*t+phi)<8*pi
    i=i+1;
    t=t+dt0;
    xb_old=xb;
    yb_old=yb;
    modom=mod(floor((om*t+phi)/pi),2);
    arg(i)=om*t;
    valmod(i)=modom;
    if modom==0
        if modom_old==1
            % ANY TIME modom SWITCHES FROM 1 TO 0 THE PHASE IS REDEFINED
            phi1=om2*(t-dt0)+phi2-om1*(t-dt0);
        end
        [xb,yb,ub,vb,theta,dtheta]=moving_coord2D_phi(n,R,t,om1,phi1);
        % MINIMUM DISTANCE BETWEEN THE PRESENT AND PREVIOUS CONFIGURATION
        % OF THE POINTS ON THE BOUNDARY, THIS IS TO MONITOR THE SPEED
        % CHANGE SINCE dt IS KEPT CONSTANT
        distb(i)=min(sqrt((xb-xb_old).^2+(yb-yb_old).^2));
        om=om1;
        phi=phi1;
    else
        if modom_old==0
            % ANY TIME modom SWITCHES FROM 0 TO 1 THE PHASE IS REDEFINED
            phi2=om1*(t-dt0)+phi1-om2*(t-dt0);
        end
        [xb,yb,ub,vb,theta,dtheta]=moving_coord2D_phi(n,R,t,om2,phi2);
        distb(i)=min(sqrt((xb-xb_old).^2+(yb-yb_old).^2));
        om=om2;
        phi=phi2;
    end
    modom_old=modom;
    time(i)=t;
    
    figure(1)
    plot(xb, yb, 'k-.')
    title(['t= ' num2str(t)])
    axis([-4 4 -4 4])
    
end

figure
plot(time,arg,'.')
hold on
plot(time,valmod,'r.')

figure
plot(time,distb,'.-')
